Molecular Simulations of Shock to Detonation Transition in Nitromethane 
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An extension of the model described in a previous work of Maillet, Soulard and Stoltz [1] based 
on a Dissipative Particule Dynamics is presented and applied to liquid nitromethane. Large scale 
non-equilibrium simulations of reacting nitromethane under sustained shock conditions allow a bet- 
ter understanding of the shock-to-detonation transition in homogeneous explosives. Moreover, the 
propagation of the reactive wave appears discontinuous since ignition points in the shocked material 
can be activated by the compressive waves emitted from the onset of chemical reactions. 

PACS numbers: 



Since two decades, molecular dynamics (MD) has be- 
come a widely used numerical tool to model and under- 
stand shock waves processes. It has been successfully 
used to study the appearance of plasticity and phase 
transition under shock conditions [H [3], as well as the 
structure of the shock front [3]. Similar successes have 
not been obtained in the study of high energetic mate- 
rials, in particular for the formation of a reactive (deto- 
nating) wave. This is mainly due to the very large time 
and length scales required to observe such phenomena. 
Simulations of reactive waves have been limited so far 
to model materials, mainly using the REBO potential or 
some of its variants O |6]. For real energetic materials, 
only the first steps of the decomposition mechanism of 
RDX under shock conditions were studied using reactive 
classical MD j7j . The same model allowed the simulation 
of the onset of the reactive wave [S]. Finally, modified 
equations of motion [i.e. constrained dynamics) allowed 
to model shock properties of reacting systems [9j , but the 
behavior of the wave itself is not accessible. In conclu- 
sion, full atomistic simulations of detonations in realis- 
tic explosives have never been achieved so far since the 
involved computational burden outperforms all current 
computer resources by orders of magnitude. 

Having in mind that detonation ultimately is a multi- 
scale process coupling hydrodynamical features and local 
chemical reactions, coarse graining strategies appear as a 
relevant alternative to a full direct numerical simulation 
of the entire detonation process. The coarse-graining ap- 
proach was only recently considered for shock waves |10l - 
[T^ . It is based on the replacement of complex molecules 
by mesoparticules with internal degrees of freedom. 

In this letter, we present equations of motions for 
mesoparticles including energy exchange between in- 
tramolecular and intermolecular degrees of freedom in 
order to ensure energy conservation, DPDE (standing 
for 'Dissipative Particle Dynamics at constant Energy'). 
Moreover, additional variables are added to each particle 



to account for its possible chemical evolution. The dy- 
namics of these variables is governed by classical kinetic 
equations. The final model, RDPDE (standing for 'Re- 
active DPDE'), is used to study the shock to detonation 
transition mechanism in nitromethane. 

Description of the model. We first briefly recall the 
DPDE model used to simulate inert shock waves in J12j . 
For a system of N particles q^ with momenta Pi = rriiVi, 
the equations of motion read 

dqi = — dt, 
rrii 

dpi = - ^q,V{q) dt - ^ijX^{r,j)vij dt + axir^j) dW,j, 

I ^^ ^, . f o 30-2/1 1 

2 \ nij rrii 



d£^ =^ 2^ X {r,j) [lijVij 
- ax{rij)vij ■ dWij, 



dt 



(1) 



where V{q) is the interaction potential between mesopar- 



ticules (a sum of paiwise interactions here). 



-m 



and Vj 



are respectively the distance and 



the relative velocity between particles i and j, x(^) is 
a weighting function with cut-off radius rcut: a-nd the 
processes Wij are independent d-dimensional Brownian 
motions with Wij — —Wji. The following fluctuation- 
dissipation relating the magnitude of the stochastic forces 
a to that of the friction forces ^ij > ensures that the 



canonical distribution is invariant: "fij 



cr^fl.: 
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2^ Ijr + ijrj. In the latter equation, the tempera- 
tures Ti are internal temperatures obtained from an in- 
ternal equation of state (EOS) on the internal energy: 

^i ~ lo' ^'"i (-^) '^"^' "^^^ explicit temperature depen- 
dence of the function C^,. allows to accurately reproduce 
the thermodynamic behavior of the internal degrees of 
freedom (DoFs). As discussed in [13], this function is 
known to modify the response of materials to shock load- 
ing. For reactive material, or in the more general case 



when some processes are activated by the temperature 
of internal DoFs, this function is expected to play an 
important role as well. 

For reactive materials, additional variables per particle 
Aj are introduced to model the progress of the chemical 
reactions at hand [1] . Their time evolution equations de- 
pend on the order of the chemical reactions and the re- 
versible natures of the latter. We have chosen to model 
the chemical decomposition of nitromethane with a suc- 
cession of two first order chemical reactions: The first one 
is reversible and endothermic, while the second one is ir- 
reversible and exothermic. One mesoparticle could repre- 
sent successively all states beteween NiME and products. 
The model reaction of the decomposition of nitromethane 



NiMe ^ NiMe* -> Products 

where NiMe* represents a metastable material corre- 
sponding to a local minimum on the potential energy 
surface. The progress variables Ai and A2 represent the 
evolutions of the two chemical reactions, and can be iden- 
tified with the fractions of NiMe* and Products respec- 
tively. Their time evolution reads: 



dXi 

~dt 

dX2 

dt 



=ki (1 - Ai - A2) - fc_iAi - fcsAi, (2) 

=feAi, (3) 



where fci and fc_i are the forward and backward reaction 
rates associated to the first chemical reaction and fc2 is 
the reaction rate associated to the second chemical reac- 
tion. The expressions of these rates are given by standard 
Arrhenius expressions : k = Aexp{—Ea/RT) where Ea 
is the activation energy and T is obtained from a local 
spacial average of internal temperatures. 

The model is constructed so that the total energy of 
the system (the sum of the potential, kinetic, internal and 
chemical energies) is constant. In particular, the energy 
variations due to chemical reactions are compensated by 
appropriate variations of the kinetic and internal ener- 
gies, as described in [T]. 

Application to nitromethane. The interactions be- 
tween mesoparticules are described by a classical exp-6 
potential, which is known to ensure a good compress- 
ibility of the system at high pressure. By means of the 
analytical expression of the equation of state of exp-6 flu- 
ids proposed by Kataoka [M] , the parameters of the po- 
tential for inert nitromethane and its detonation prod- 
ucts are obtained respectively by optimization on ex- 
perimental Hugoniot [TMT7] and Crussard curves [TH]. 
Consistency of these potentials with experimental data 
is shown on Figure [T] From now on, computations are 
made with these potentials. When a DPD particule is 
reacting (0 < Aj < 1), the parameters of the interaction 
potential are given by a linear interpolation of the param- 



eters of the potential for the inert and the fully reacted 
system. 

For inert nitromethane, as the DPD particle represents 
only a single nitromethane molecule, there is a clear split- 
ting between internal and external DoFs: the internal 
Cv of the DPD particle represents simply the sum of the 
intramolecular contributions to the energy. Hence the 
Cy function is taken directly from thermodynamic ta- 
bles. For the fully reacted material, a DPD particle rep- 
resents the decomposition products of one nitromethane 
molecule, i.e. a group of several small molecules. The 
corresponding internal C„ is then a sum of two contri- 
butions: the internal DoFs of the small molecules plus a 
fraction of their external DoFs. There is no clear split- 
ting between intra and intermolecular DoFs in the coarse 
graining process. Therefore, the C^ function of the DPD 
particle cannot be extracted from a standard thermody- 
namic table, and has to be computed numerically as the 
difference between the total Cy of the real system and 
the Cv due to the intermolecular DPD potential only. 
In the following internal DoFs will stands for all DoFs 
inside a mesoparticle and external DoFs will stands for 
displacement of mesoparticles. 

The rate of heat exchanges between internal and ex- 
ternal DoFs is controlled by the parameter a. Dawes 
et al [19] have presented a study of energy relaxation in 
solid nitromethane behind a shock wave. They computed 
the rate of energy transfer between translational and ro- 
tational kinetic energies and internal vibrational modes. 
Of course, some of the vibrational modes that exhibit 
an efficient coupling to the phonons are therefore heated 
rapidly after the shock, while some other modes heat up 
more slowly. In the DPD model we use, there is only one 
mode representing the average of all internal vibrations. 
We chose a value of a leading to an equilibration time of a 
few picoseconds between intra and intermolecular DoFs. 

Numerical results. All simulations have been per- 
formed with our in-house parallel MD code. We used 
a time step At = 1 fs in our simulations, a value at least 
one order of magnitude larger than for MD simulations 
with standard reactive potentials as such ReaxFF. 

The thermodynamic properties of our model have been 
studied by equilibrium and nonequilibrium MD simula- 
tions. The computed density of the material at P = 1 bar 
and T — 300 K is 1.1041g.cm^^. Simulations of inert ni- 
tromethane using the hugoniostat constraint method and 
NEMD (shock simulations) have been performed in or- 
der to compute the inert Hugoniot curve. Similar simu- 
lations are performed with reacted nitromethane to com- 
pute the Hugoniot curve of detonation products, i.e. the 
Crussard curve. Results are displayed in Figure [T] A 
semi-quantitative agreement is found on thermodynami- 
cal properties of both the inert and reacted explosive. 

Nonequilibrium shock wave simulations have been per- 
formed with the RDPDE model. Sustained shock simu- 
lations have been loaded on sample of size 20x20x 10000 
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FIG. 1: Hugoniot and Crussard curves of nitromethane. 
Black line: experimental Hugoniot curve. Dotted line: Crus- 
sard curve (CARTE code). Open symbols correspond to in- 
ert nitromethane while full symbols correspond to detona- 
tion products. Triangles: calculations made with the EOS 
of Kataoka. Squares: Monte Carlo simulations. Diamonds: 
MD-DPD equilibrium simulations. Circles: non equilibrium 
shock simulations. 



unit cells (about 4.5/im long), with an infinitely massive 
piston moving at constant velocities equal to 1500, 2000, 
2500, 2700 and 3000 m.s'^ The total simulation time 
for each simulation is around 700 ps. In the simulations 
with piston velocities equal to 1500 and 2000 m.~^, no 
sign of reactive wave formation is observed. The inert 
shock wave propagates through the whole sample, acti- 
vating the explosive with the first reversible reaction, but 
the detonation products are not produced. This is due to 
the fact that the shock did not release enough energy for 
the system to overcome the energy barrier of the second 
(irreversible) reaction. For piston velocities equal to or 
higher than 2500 m.s~^, we observe the decomposition of 
the explosive and the formation of a reactive wave. The 
velocity of the corresponding wave now exceeds the one 
of the inert shock wave due to the global energy release 
of exothermic reactions. 

On the experimental side, multiple-magnetic gauge 
measurements of neat nitromethane [2D] allowed the di- 
rect observation of the shock-to-detonation mechanism: 
A reactive wave first builds up in the shocked (inert) ma- 
terial (characterized by an increase of the particle veloc- 
ity as the wave propagates) and forms a superdetonation 
wave travelling at constant velocity. This superdetona- 
tion wave overtakes the inert shock wave and progres- 
sively decays to an overdriven, and eventually to a steady, 
detonation. This has recently been confirmed by PDV 
(Photonic Doppler Velocimetry) measurements [21j . 

For the simulation using a piston velocity of 
2700 m.s^^, the onset of the reactive wave is evi- 
denced in Figure [2] by the signal of (artificial) multi- 
ple gauges. This wave progressively catches up with the 
shock wave. These signals exhibit striking similarities 




FIG. 2: Time evolutions of particle velocities taken at differ- 
ent locations in the material (every 50 cells), mimicking the 
behaviour of multiple gauges set up. 



with experimental measurements, although occuring on 
a shorter timescale. The complete scenario of the shock- 
to-detonation transition is displayed in Figure |3l which 
shows the evolutions of the average velocity, progress 
variable, internal temperature and the x-component of 
the pressure tensor, in a (x, t) diagram. 

After a shock wave is loaded in the sample it leaves the 
nitromethane in an shocked state. After around 100 ps, 
chemical reactions start to occur at the interface with the 
piston since this is the region of space where the system 
has spent the longest time in the shocked state. As the 
chemical reactions progress, a reactive wave is formed, 
which travels in the shocked nitromethane at a higher 
velocity than the shock wave. Its front is given by the 
analysis of the diagram of the progress variable. From the 
onset of chemical reactions happening around 100 ps, a 
compressive wave is created, travelling much faster than 
the reactive wave in the shocked material, up to the in- 
ert shock front. This compressive wave brings the inert 
shocked system at a slightly higher density. When it 
reaches the shock front, it accelerates and increases the 
pressure (and the material velocity) in the shocked in- 
ert material. At some point, an ignition zone induced 
by the previously described compressive wave appears in 
the region between the inert shock wave and the reac- 
tive wave. This ignition point acts as a hot spot, and 
chemical reactions develop from this point forward and 
backward in the material, forming two reacting waves. 
Once again, these reacting waves are preceded by com- 
pressive waves. A careful analysis of the diagram reveals 
that a second hot spot is activated, closer to the front 
shock. As the reactive wave develops from this point, 
it catches up with the shock wave, forming a single re- 
active wave, the overdriven detonation. This overdriven 



4»4> 


^^ 




asSE- 


-1»M 


^^^^M 




M» 


rl^S^Ofl 


^^^^^H 




M0& 


■3*-oa 


^"^^^^^^^^^^Hl 






-£4-^7 


1 j^B 




y.ij 


-it-se 


^hhI 




afi 


•l.5t<* 


^^^^^^^^^^^^H 


■ 




-3*-09 


Hi^^^^^l 


1 


0? 




~^^^^^^^^^B 


.. 


tritiL- 


^SUT 


^<^flfl 




*«* 


-it4ft 


^^S 




3SW 


■1AH» 


.^fl^^l 




£006 


*« 


■^^^^^^^^^^H 




IME- 


4*<fr 


^2 






-1*-M 


jt^^^k 




£Se*1d 


-1,ta-» 


.^fl^^l 




I.Sa*1Q 


-&»-W 


;^^^^^^^^^H 




1*+ICl 



discussions. G.S. is supported in part by the Agence Na- 
tionale de la Recherche, under the grant ANR-10-BLAN 
0108 (MEG AS). All Monte Carlo simulations have been 
performed with the GIBBS code from IFF, CNRS and 
the Universite Faris-Sud 122 . 
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FIG. 3: Coloured maps of thermodynamic variable (particle 
velocity, progress variable, internal temperature and pressure) 
in a time-space diagram. 



detonation propagates at high speed and is characterized 
by high particle velocities, as well as high values of ther- 
modynamic variables (temperature, pressure,...). Once 
formed, the overdriven detonation progressively decays 
to a stationary detonation, as evidenced in the pressure 
and temperature diagrams. 

In conclusion, we presented for the first time a molec- 
ular simulation of a shock-to-detonation transition in a 
realistic energetic material, using a reduced model based 
on a DFD coarse graining approach. Our numerical re- 
sults confirmed the shock to detonation mechanism ob- 
served experimentally for nitromethane. Moreover, we 
have shown that compressive waves are emitted from the 
onset of chemical reactions, therefore increasing the prob- 
ability of ignition forward in the material. Hence the 
reactive wave exhibits a discontinuous progression from 
sites to sites up to the shock front. 

Laurent Soulard is gratefully acknowledged for fruitful 
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